## REPLICATION FOR: "From Garments to Grievances: The Dynamics of Urban Growth and Economic Shocks in Rapidly Industrializing Economies"
## Author: Elisa D'Amico
## Date: May 19, 2025


rm(list=ls())

library(data.table)
library(did)
library(dplyr)
library(fixest)
library(ggplot2)
library(haven)
library(knitr)
library(tidyverse)

# Import Data
# -----------
H1DID <- read.csv("DID.csv")
H3DDD <- read.csv("DID.csv")

# MAIN ANALYSIS: HYPOTHESIS 1
# ---------------------------

# Main model
# -------------
H1result_main <- att_gt(
  yname = "urbanproportion",
  tname = "year",
  idname = "admnum",
  gname = "firsttreat",
  data = H1DID,
  bstrap = TRUE,
  control_group = "notyettreated"
)

# Dynamic event study estimates
csdid_stats_event_main <- aggte(H1result_main, type = "dynamic", na.rm = TRUE)
summary(csdid_stats_event_main)

# Plot main results
main_plot <- ggdid(csdid_stats_event_main) +
  labs(
    # title = "DID Event Study: Economic Demand Shocks and Urban Migration",
    title = "",
    x = "Time to Treatment (Years)",
    y = "Average Treatment Effect on Urban Proportion"
  ) +
  scale_color_manual(
    values = c("#9B59B6", "#00BFC4"),
    labels = c("Pre-Treatment", "Post-Treatment")
  ) +
  theme(
    plot.title = element_text(size = 14, face = "bold", hjust = 0.5, family = "Times New Roman"),
    axis.text = element_text(size = 8, family = "Times New Roman"),
    axis.title = element_text(size = 10, family = "Times New Roman"),
    legend.text = element_text(family = "Times New Roman"),
    legend.title = element_text(family = "Times New Roman")
  )

print(main_plot)

# ROBUSTNESS CHECKS
# -----------------

# 1. Parallel Trends Assumption Check
# -----------------------------------
coefs <- c(0.0000, 0.0000, 0.0005, 0.0000, 0.0000, 0.0002, 0.0008, 0.0009, -0.0004,
           0.0004, 0.0000, -0.0002, 0.0027, 0.0009, 0.0035, 0.0026, 0.0024, 0.0094,
           -0.0060, -0.0052, 0.0102, 0.0065, 0.0024)

lower_ci <- c(0.0000, 0.0000, -0.0003, 0.0000, 0.0000, -0.0005, -0.0010, -0.0020, -0.0019,
              -0.0014, -0.0013, -0.0025, -0.0002, -0.0009, -0.0008, -0.0005, -0.0042, -0.0054,
              -0.0174, -0.0190, -0.0057, -0.0067, -0.0082)

upper_ci <- c(0.0000, 0.0000, 0.0014, 0.0000, 0.0000, 0.0008, 0.0026, 0.0037, 0.0011,
              0.0023, 0.0014, 0.0020, 0.0057, 0.0027, 0.0079, 0.0057, 0.0089, 0.0242,
              0.0054, 0.0085, 0.0261, 0.0197, 0.0131)

time_periods <- -22:0


# Figure D1. Pretrend coefficients in urbanization between RMG factory emergent and control group cities (95% CI)
# Create pre-trends plot
plot(time_periods, coefs, type = "n", ylim = c(min(lower_ci), max(upper_ci)),
     xlab = "Time Period", ylab = "Coefficient",
     main = "")
     # main = "Pre-Treatment Coefficients (95% CI)")

segments(time_periods, lower_ci, time_periods, upper_ci, col = "black", lwd = 2)
abline(h = 0, col = "red", lty = 2)
points(time_periods, coefs, pch = 19, col = "black")


# 2. Inverse Probability Weighting
# --------------------------------
H1result_ipw <- att_gt(
  yname = "urbanproportion",
  tname = "year",
  idname = "admnum",
  gname = "firsttreat",
  data = H1DID,
  est_method = "ipw"
)

csdid_stats_event_ipw <- aggte(H1result_ipw, type = "dynamic", na.rm = TRUE)
summary(csdid_stats_event_ipw)
print(ggdid(csdid_stats_event_ipw) +
        labs(title = "DID Event Study: Economic Shocks and Urban Migration (IPW)"))

# 3. Clustered Standard Errors
# ----------------------------
H1result_cluster <- att_gt(
  yname = "urbanproportion",
  tname = "year",
  idname = "admnum",
  gname = "firsttreat",
  data = H1DID,
  clustervars = c("admnum")
)

csdid_stats_event_cluster <- aggte(H1result_cluster, type = "dynamic", na.rm = TRUE)
summary(csdid_stats_event_cluster)
print(ggdid(csdid_stats_event_cluster) +
        labs(title = "DID Event Study: Economic Shocks and Urban Migration (Clustered SEs)"))

# 4. Custom Alpha (0.01)
# ----------------------
H1result_custom_alpha <- att_gt(
  yname = "urbanproportion",
  tname = "year",
  idname = "admnum",
  gname = "firsttreat",
  data = H1DID,
  alp = 0.01
)

csdid_stats_event_custom_alpha <- aggte(H1result_custom_alpha, type = "dynamic", na.rm = TRUE)
summary(csdid_stats_event_custom_alpha)
print(ggdid(csdid_stats_event_custom_alpha) +
        labs(title = "DID Event Study: Economic Shocks and Urban Migration (α=0.01)"))

# 5. Never Treated Control Group
# ------------------------------
H1result_never_treated <- att_gt(
  yname = "urbanproportion",
  tname = "year",
  idname = "admnum",
  gname = "firsttreat",
  data = H1DID,
  control_group = "nevertreated"
)

csdid_stats_event_never_treated <- aggte(H1result_never_treated, type = "dynamic", na.rm = TRUE)
summary(csdid_stats_event_never_treated)
print(ggdid(csdid_stats_event_never_treated) +
        labs(title = "DID Event Study: Economic Shocks and Urban Migration (Never Treated Controls)"))

# 6. Not Yet Treated Control Group (Redundant but included)
# --------------------------------------------------------
H1result_not_yet_treated <- att_gt(
  yname = "urbanproportion",
  tname = "year",
  idname = "admnum",
  gname = "firsttreat",
  data = H1DID,
  control_group = "notyettreated"
)

csdid_stats_event_not_yet_treated <- aggte(H1result_not_yet_treated, type = "dynamic", na.rm = TRUE)
summary(csdid_stats_event_not_yet_treated)
print(ggdid(csdid_stats_event_not_yet_treated) +
        labs(title = "DID Event Study: Economic Shocks and Urban Migration (Not Yet Treated Controls)"))


# MAIN ANALYSIS: HYPOTHESIS 2
# ---------------------------

# Figure 6. Impact of RMG factory presence on urban unrest events.
# Direct Migration Unrest Relationship

# Estimate DID model
H2result_DirectMigUnrest <- att_gt(
  yname = "urbanunrestcount",
  tname = "year",
  idname = "admnum",
  gname = "firsttreat",
  data = H1DID,
  bstrap = TRUE,
  control_group = "notyettreated"
)

csdid_stats_event_migun <- aggte(H2result_DirectMigUnrest, type = "dynamic", na.rm = TRUE)
summary(csdid_stats_event_migun)

# Plot event study results
unrest_plot <- ggdid(csdid_stats_event_migun) +
  labs(
    title = "",
    x = "Time to Treatment (Years)",
    y = "Average Treatment Effect (ATT)"
  ) +
  scale_color_manual(
    values = c("#6A0DAD", "#00CED1"),
    labels = c("Pre-Treatment", "Post-Treatment")
  ) +
  theme(
    plot.title = element_text(size = 14, face = "bold", hjust = 0.5, family = "Times New Roman"),
    axis.text = element_text(size = 8, family = "Times New Roman", color = "black"),
    axis.title = element_text(size = 10, family = "Times New Roman", color = "black"),
    legend.text = element_text(family = "Times New Roman"),
    legend.title = element_text(family = "Times New Roman"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

print(unrest_plot)


# MAIN ANALYSIS: HYPOTHESIS 3
# ---------------------------

# Prepare data for DDD analysis
# -----------------------------
subset_H3DDD <- H3DDD %>%
  filter(firsttreat != 0) %>%
  arrange(admnum, year) %>%
  group_by(admnum) %>%
  slice(1) %>%
  ungroup() %>%
  mutate(firstfact_binary = 1) %>%
  select(admnum, year, firstfact_binary)

H3DDD <- H3DDD %>%
  left_join(subset_H3DDD, by = c("admnum", "year")) %>%
  mutate(firstfact_binary = replace_na(firstfact_binary, 0))

setDT(H3DDD)
H3DDD[, treat := as.integer(!is.na(firsttreat) & firsttreat != 0)]
H3DDD[, time_to_treat := fifelse(treat == 1, year - firsttreat, 0)]
H3DDD[, never_treated := as.integer(firsttreat == 0)]


H3DDD <- H3DDD %>%
  arrange(admnum, year) %>%
  group_by(admnum) %>%
  mutate(
    log_foodprices_global = log(foodprices_global),
    interaction_term = firstfact_binary * log_foodprices_global,
    urbanunrest_binary = as.integer(urbanunrestcount > 0),
    urbanunrest_binary_lag = lag(urbanunrest_binary)
  ) %>%
  ungroup()

# MAIN MODEL: H3
# ---------------

# Subset data around the treatment period
H3DDDsub <- H3DDD %>%
  filter(time_to_treat >= -3 & time_to_treat <= 3)

# Run the TWFE event study model
mod_twfe <- feols(urbanunrest_binary + never_treated + interaction_term + urbanproportion ~ i(time_to_treat, treat, ref = -1) +
                    interaction_term + never_treated + urbanproportion |
                    year + admnum,
                  cluster = ~admnum + never_treated,
                  data = H3DDDsub)

summary(mod_twfe)

# Extract ATT coefficients
att_coefficients <- c(
  coef(mod_twfe)["time_to_treat::1:treat"],
  coef(mod_twfe)["time_to_treat::2:treat"],
  coef(mod_twfe)["time_to_treat::3:treat"]
)

average_ATT <- mean(att_coefficients)
average_ATT

# Figure 7. DDD Event plot: negative economic shocks X urbanization on unrest.
# Plot TWFE event study results
iplot(mod_twfe,
      xlab = 'Time to Treatment',
      ylab = 'Average Treatment Effect (ATT)',
      # main = 'DDD Event Plot: Negative Economic Shocks X Migration and Unrest')
      main = '')
abline(v = -1, lty = 2, col = "white")
abline(v = 0, lty = 2, col = "black")

# SUN AND ABRAHAM (Robustness Check)
# ----------------------------------

# Create year_treated variable for Sun & Abraham method
setDT(H3DDDsub)
H3DDDsub[, year_treated := ifelse(treat == 0, 10000, `firsttreat`)]

# Run Sun and Abraham event study model
mod_sa = feols(urbanunrest_binary + never_treated + interaction_term + urbanproportion ~ sunab(year_treated, year) +
                 interaction_term + never_treated + urbanproportion |
                 year + admnum,
               cluster = ~admnum + never_treated,
               data = H3DDDsub)

summary(mod_sa)

# Figure N1. EPlotted Figure 7 with TWFE Robustness.
# Compare TWFE and Sun & Abraham results
iplot(list(mod_twfe, mod_sa), sep = 0.5, ref.line = -1,
      xlab = 'Time to Treatment',
      ylab = 'Average Treatment Effect (ATT)',
      main = '')
    #  main = 'DDD Event Plot: Negative Economic Shocks X Migration on Unrest')
legend("bottomleft", col = c(1, 2), pch = c(20, 17),
       legend = c("Sun & Abraham (2020)", "TWFE"))
abline(v = -1, lty = 2, col = "white")
abline(v = 0, lty = 2, col = "black")


# Parallel Trends Plot (Figure F1)
# --------------------------------

# Data for the plot
time_periods <- c(-20:-19,-17:-1,0)
coefficients <- c(-0.09897, -0.093887, -0.090259, -0.08094, -0.079226, -0.075882, 
                  -0.075742, -0.071233, -0.071804, -0.070372, -0.066997, -0.061572, 
                  -0.059439, -0.05684, -0.053371, -0.055888, -0.033583, 0.001464, 
                  -0.008558, 0.007677)
lower_ci <- c(-11.338543, -10.892774, -10.048236, -9.463509, -9.719178, -10.067291, 
              -10.611359, -10.129188, -10.476518, -11.448155, -13.251323, -14.023589, 
              -23.432629, -32.40316, -49.820162, -1130.704337, -20.572304, -0.741665, 
              -3.38589, -10.569485)
upper_ci <- c(0.056002, 0.058281, 0.063148, 0.067022, 0.065272, 0.063, 0.059818, 
              0.062647, 0.060583, 0.055468, 0.047951, 0.04532, 0.027152, 0.019641, 
              0.012777, 0.000563, 0.030921, 0.59374, 0.18282, 0.060053)

y_min <- -0.1
plot(time_periods, coefficients, type = "n", ylim = range(y_min, upper_ci),
     xlab = "Time Period", ylab = "Coefficient", 
     main = "")
     # main = "Pretreatment Coefficients (95% CI)")
segments(time_periods, lower_ci, time_periods, upper_ci, col = "black", lwd = 2)
abline(h = 0, col = "red", lty = 2)
points(time_periods, coefficients, pch = 19, col = "black")


# Figure F1. Pretrend coefficients in urban unrest between RMG factory emergent and control group cities contingent on food prices (95% CIs)
coef_table <- data.frame(Time_Period = time_periods,
                         Coefficient = coefficients,
                         Lower_CI = lower_ci,
                         Upper_CI = upper_ci)
print(coef_table)

